mydata <- read.csv2("People analytics.csv")
head(mydata)
##   ID             Track Econ Comm Analyt Finance
## 1  1 Operational track  4.6  4.8    3.5     3.1
## 2  2 Operational track  3.7  3.4    2.7     3.1
## 3  3 Operational track  4.4  4.8    3.3     3.3
## 4  4 Operational track  4.5  4.7    2.5     3.4
## 5  5 Operational track  3.9  3.5    3.4     3.3
## 6  6 Operational track  3.7  3.0    3.4     4.0

Description of the variables:

mydata$Track <- factor(mydata$Track)
mydata_DA <- mydata[c(3:6)] 

desc_stats <- data.frame(
  N        = apply(mydata_DA, 2, function(x) sum(!is.na(x))),
  Missing  = apply(mydata_DA, 2, function(x) sum(is.na(x))),
  Mean     = apply(mydata_DA, 2, mean, na.rm = TRUE),
  Variance = apply(mydata_DA, 2, var, na.rm = TRUE),
  Min      = apply(mydata_DA, 2, min, na.rm = TRUE),
  Max      = apply(mydata_DA, 2, max, na.rm = TRUE)
)

round(desc_stats, 2)
##           N Missing Mean Variance Min Max
## Econ    310       0 4.29     0.48 2.0 6.0
## Comm    310       0 4.26     0.55 2.3 6.0
## Analyt  310       0 4.23     0.62 2.4 6.0
## Finance 310       0 4.03     0.54 2.2 5.6
table(mydata$Track)
## 
##  Analytical track Operational track   Strategic track 
##               105                95               110
# Multiple ANOVAs by Track
fit <- aov(cbind(Econ, Comm, Analyt, Finance) ~ Track, data = mydata)
summary(fit)
##  Response Econ :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Track         2  11.07  5.5351  12.497 6.058e-06 ***
## Residuals   307 135.98  0.4429                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Comm :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Track         2 76.494  38.247  127.34 < 2.2e-16 ***
## Residuals   307 92.208   0.300                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Analyt :
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## Track         2  86.461  43.230  125.19 < 2.2e-16 ***
## Residuals   307 106.017   0.345                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Finance :
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## Track         2  55.559 27.7797  77.554 < 2.2e-16 ***
## Residuals   307 109.966  0.3582                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Multicolinearity
vars <- c("Econ", "Comm", "Analyt", "Finance")

vif_all(mydata, vars)
## $Econ
##     Comm   Analyt  Finance 
## 1.001786 1.274940 1.274983 
## 
## $Comm
##     Econ   Analyt  Finance 
## 1.061216 1.314375 1.280724 
## 
## $Analyt
##     Econ     Comm  Finance 
## 1.032028 1.004366 1.030999 
## 
## $Finance
##     Econ     Comm   Analyt 
## 1.059498 1.004668 1.058407
# Box’s M test for homogeneity of covariance matrices across Track groups
library(biotools)
boxM(mydata[, c("Econ", "Comm", "Analyt", "Finance")],
     mydata$Track)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  mydata[, c("Econ", "Comm", "Analyt", "Finance")]
## Chi-Sq (approx.) = 22.064, df = 20, p-value = 0.3371
# Computes squared Mahalanobis distance to identify multivariate outliers based on variable covariance structure.
mydata <- mahal_distance(
  data = mydata,
  vars = c("Econ", "Comm", "Analyt", "Finance"),
  group = "Track",
  prob = 0.99)

outliers <- mydata %>% 
  filter(outlier)

outliers
## # A tibble: 3 × 9
##   Track             ID  Econ  Comm Analyt Finance    md cutoff outlier
##   <fct>          <int> <dbl> <dbl>  <dbl>   <dbl> <dbl>  <dbl> <lgl>  
## 1 Analytical tr…   199   3     5.1    5.8     3.4  13.7   13.3 TRUE   
## 2 Operational t…    41   5.5   4.7    5.2     4.4  13.5   13.3 TRUE   
## 3 Strategic tra…   290   5.6   3.8    4.4     2.3  15.2   13.3 TRUE
# Creates a Chi-square Q-Q plot to compare observed squared Mahalanobis distances with theoretical values for assessing multivariate normality.
n <- length(mydata$md)
theoretical <- qchisq(ppoints(n), df = 4)
observed <- sort(mydata$md)

plot(theoretical, observed,
     main = "Chi-square Q-Q plot (Squared Mahalanobis distance)",
     xlab = "Theoretical Chi-square quantiles",
     ylab = "Observed squared Mahalanobis distance",
     pch = 20)
abline(0, 1, col = "darkgreen", lwd = 2)

# Creates separate Chi-square Q-Q plots to assess multivariate normality within each Track group
md_by_group <- split(mydata$md, mydata$Track)

par(mfrow = c(1, 3))

for (g in names(md_by_group)) {
  
  observed <- sort(md_by_group[[g]])
  n <- length(observed)
  theoretical <- qchisq(ppoints(n), df = 4)
  
  plot(
    theoretical, observed,
    main = g,
    xlab = "Theoretical Chi-square quantiles",
    ylab = "Observed squared Mahalanobis distance",
    pch = 20,
    asp = 1
  )
  
  abline(0, 1, col = "darkgreen", lwd = 2)
}

par(mfrow = c(1, 1))
#Linear discriminannt analysis
lda_model <- lda(Track ~ Econ + Comm + Analyt + Finance,
                 data = mydata)
lda_model
## Call:
## lda(Track ~ Econ + Comm + Analyt + Finance, data = mydata)
## 
## Prior probabilities of groups:
##  Analytical track Operational track   Strategic track 
##         0.3387097         0.3064516         0.3548387 
## 
## Group means:
##                       Econ     Comm   Analyt  Finance
## Analytical track  4.027619 3.994286 4.908571 4.562857
## Operational track 4.466316 3.789474 3.606316 3.514737
## Strategic track   4.380000 4.920909 4.127273 3.960909
## 
## Coefficients of linear discriminants:
##                LD1        LD2
## Econ    -0.2880613 0.04679253
## Comm    -0.2016790 1.81157732
## Analyt   1.2518192 0.05949666
## Finance  0.9100927 0.10336903
## 
## Proportion of trace:
##    LD1    LD2 
## 0.5963 0.4037
# Extract eigenvalues to assess the importance of LDA discriminant functions
eigenvalues <- lda_eigenvalues(
  mydata   = mydata,
  group    = "Track",
  variables = c("Econ", "Comm", "Analyt", "Finance")
)

eigenvalues
## [1] 1.2264305 0.8303744
#Canonical correlations and squared canonical correlations
round(canonical_corr_from_eigen(eigenvalues), 3)
##     CanonicalCorr SqCanonicalCorr
## LD1         0.742           0.551
## LD2         0.674           0.454
# Multivariate ANOVA (MANOVA) testing Track effects on all outcomes jointly
man <- manova(cbind(Econ, Comm, Analyt, Finance) ~ Track,
              data = mydata)
# Canonical discriminant analysis based on the MANOVA model
library(candisc)
can <- candisc(man)
can
## 
## Canonical Discriminant Analysis for Track:
## 
##    CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.55085    1.22643    0.39606  59.628     59.628
## 2 0.45366    0.83037    0.39606  40.372    100.000
## 
## Test of H0: The canonical correlations in the 
## current row and all that follow are zero
## 
##   LR test stat approx F numDF denDF   Pr(> F)    
## 1      0.24539   77.422     8   608 < 2.2e-16 ***
## 2      0.54634   84.421     3   305 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Raw (unstandardized) canonical discriminant coefficients (weights)
can$coeffs.raw
##               Can1       Can2
## Econ     0.2880613 0.04679253
## Comm     0.2016790 1.81157732
## Analyt  -1.2518192 0.05949666
## Finance -0.9100927 0.10336903
# Standardized canonical discriminant coefficients (beta weights)
can$coeffs.std
##               Can1       Can2
## Econ     0.1917124 0.03114166
## Comm     0.1105289 0.99282312
## Analyt  -0.7356302 0.03496315
## Finance -0.5446863 0.06186589
# Structure matrix (canonical loadings)
DA_structure_matrix(mydata   = mydata,
                    group    = "Track",
                    lda_model = lda_model,
                    variables = c("Econ", "Comm", "Analyt", "Finance"))
##                 LD1        LD2
## Econ    -0.25534686 0.04173591
## Comm    -0.06028717 0.99683350
## Analyt   0.81204226 0.09059543
## Finance  0.63710169 0.09461488
# Computes the discriminant scores (LD1, LD2, ...) by projecting each observation onto the linear discriminant functions
lda_scores <- predict(lda_model)$x

# Compute group centroids (mean LDA scores) for each Track
centroids <- data.frame(
  Track = mydata$Track,
  lda_scores) %>%
  dplyr::group_by(Track) %>%
  dplyr::summarise(across(everything(), mean))

centroids
## # A tibble: 3 × 3
##   Track                LD1    LD2
##   <fct>              <dbl>  <dbl>
## 1 Analytical track   1.46  -0.399
## 2 Operational track -1.21  -0.935
## 3 Strategic track   -0.353  1.19
# Scatterplot of LDA discriminant scores by Track
data_lda <- bind_cols(mydata, as.data.frame(lda_scores))

ggplot(data_lda, aes(LD1, LD2, color = Track, shape = Track)) +
  geom_point(size = 2.0, alpha = 0.7) +
  stat_ellipse(type = "t", level = 0.68) +
  labs(
    title = "Canonical Discriminant Functions",
    x = "Discriminant Function 1",
    y = "Discriminant Function 2"
  ) +
  theme_minimal() +
  theme(legend.position = "right")

# Confusion matrix comparing actual Track membership with LDA classifications
pred <- predict(lda_model)

table(Actual    = mydata$Track,
      Predicted = pred$class)
##                    Predicted
## Actual              Analytical track Operational track
##   Analytical track                87                 5
##   Operational track                7                79
##   Strategic track                 14                11
##                    Predicted
## Actual              Strategic track
##   Analytical track               13
##   Operational track               9
##   Strategic track                85
# LDA classification accuracy
mean(pred$class == mydata$Track)
## [1] 0.8096774
# Cross-validated LDA classification using leave-one-out validation
lda_cv <- lda(Track ~ Econ + Comm + Analyt + Finance,
              data = mydata,
              CV = TRUE)

table(Actual    = mydata$Track,
      Predicted = lda_cv$class)
##                    Predicted
## Actual              Analytical track Operational track
##   Analytical track                85                 5
##   Operational track                7                79
##   Strategic track                 15                11
##                    Predicted
## Actual              Strategic track
##   Analytical track               15
##   Operational track               9
##   Strategic track                84
# Cross-validated classification accuracy
mean(lda_cv$class == mydata$Track)
## [1] 0.8
# Cross-validated classification accuracy by group
confusin_matrix <- table(Actual = mydata$Track, Predicted = lda_cv$class)

group_accuracy <- diag(prop.table(confusin_matrix, 1))
group_accuracy
##  Analytical track Operational track   Strategic track 
##         0.8095238         0.8315789         0.7636364
# Add discriminant scores (LD1, LD2, …) and predicted group to mydata
mydata <- cbind(
  mydata,
  as.data.frame(pred$x),          # D1, D2 (LD1, LD2)
  Predicted_Track = pred$class    # predicted group membership
)
# Add posterior probabilities
mydata <- cbind(
  mydata,
  round(as.data.frame(pred$posterior), 3)
)
head(mydata, 3)
##              Track ID Econ Comm Analyt Finance       md  cutoff
## 1 Analytical track 96  3.6  3.0    6.0     4.5 8.420470 13.2767
## 2 Analytical track 97  5.0  3.9    4.6     4.0 2.995015 13.2767
## 3 Analytical track 98  5.3  4.3    4.3     4.2 4.711675 13.2767
##   outlier         LD1        LD2  Predicted_Track Analytical track
## 1   FALSE  3.09450498 -2.1613645 Analytical track            0.999
## 2   FALSE  0.30211480 -0.6004153 Analytical track            0.529
## 3   FALSE -0.05850241  0.1410782  Strategic track            0.244
##   Operational track Strategic track
## 1             0.001           0.000
## 2             0.290           0.181
## 3             0.235           0.521
new_person <- data.frame(
  Econ    = 2.0,
  Comm    = 2.7,
  Analyt  = 4.2,
  Finance = 2.9
)
pred_new <- predict(lda_model, newdata = new_person)
pred_new$class        # Assigned group
## [1] Operational track
## Levels: Analytical track Operational track Strategic track
pred_new$x            # LD1 and LD2 scores (centered and comparable)
##           LD1      LD2
## 1 -0.09351622 -3.05219
pred_new$posterior    # Posterior probabilities
##   Analytical track Operational track Strategic track
## 1        0.1454469         0.8524696     0.002083432